Local bifurcation with spin-transfer torque in superparamagnetic tunnel junctions

Modulation of the energy landscape by external perturbations governs various thermally-activated phenomena, described by the Arrhenius law. Thermal fluctuation of nanoscale magnetic tunnel junctions with spin-transfer torque (STT) shows promise for unconventional computing, whereas its rigorous representation, based on the Néel-Arrhenius law, has been controversial. In particular, the exponents for thermally-activated switching rate therein, have been inaccessible with conventional thermally-stable nanomagnets with decade-long retention time. Here we approach the Néel-Arrhenius law with STT utilising superparamagnetic tunnel junctions that have high sensitivity to external perturbations and determine the exponents through several independent measurements including homodyne-detected ferromagnetic resonance, nanosecond STT switching, and random telegraph noise. Furthermore, we show that the results are comprehensively described by a concept of local bifurcation observed in various physical systems. The findings demonstrate the capability of superparamagnetic tunnel junction as a useful tester for statistical physics as well as sophisticated engineering of probabilistic computing hardware with a rigorous mathematical foundation.

A dynamical system is classified by the stability of its potential landscape, especially by the local bifurcation of the dynamical equation. Under finite stochasticity, the Arrhenius law, a general principle for thermally-activated events, describes a variety of dynamical phenomena ranging from chemical reactions to physical processes. According to the Arrhenius law, the relaxation time τ in staying at a certain state is given by τ = τ 0 expΔ with the thermal stability factor Δ ≡ E 0 /k B T, where τ 0 is the intrinsic time constant of each system, k B the Boltzmann constant, T the absolute temperature, and E 0 an intrinsic energy barrier for switching to different states without perturbation. Under the perturbation by normalised external input x, E 0 is replaced by an effective energy barrier E = E 0 1 À x ð Þ n x , where the switching exponent n x is determined by the effective energy landscape with x. In general, when the dynamical equation dΘ/ dt = f(Θ,x) is given, where t is the time and Θ is the state variable (e.g. particle position, spin quantisation direction, amount of unreacted substance, etc.), the types of local bifurcation of f(Θ,x) determines n x ; for example, n x = 2 for the pitchfork bifurcation and n x = 3/2 for the saddle-node bifurcation 1 . This suggests that the switching exponent serves as a unique lens for the local structure, especially the stability of the energy landscape under the perturbations in the relevant system.
In magnetic materials, which have served as a model system to study the physics of thermally-activated phenomena, the basis of the Arrhenius law was built by Néel 2 and Brown 3 , known as the Néel-Arrhenius law. For single-domain uniaxial magnets with a magnetic field H applied along the easy axis, E under magnetic field can be simply derived as E = E 0 (1 -H/H K eff ) 2 by the Stoner-Wohlfarth model, where H K eff is the effective magnetic anisotropy field 4 . However, theoretical studies pointed out that the value of exponent n H , 2 in the above equation, should vary when one considers some realistic factors such as misalignment of magnetic field 5,6 and higher-order terms of anisotropy 7 ; in other words, the local bifurcation varies with them.
The magnetisation of nanomagnets can also be controlled by spin-transfer torque (STT) under current application through the angular momentum transfer [8][9][10][11][12][13] . The STT-induced magnetisation switching in thermally-stable magnetic tunnel junctions (MTJ) is a key ingredient for non-volatile magnetoresistive random access memory [14][15][16] . Moreover, recent studies have demonstrated an unconventional paradigm of computing, e.g. neuromorphic computing with population coding 17 , and probabilistic computing 18 , which utilises a combinatorial effect of STT and thermal fluctuation in superparamagnetic tunnel junctions (s-MTJs). By further combining the effect of external magnetic fields, additional tunabilities of the s-MTJs for probabilistic computing have been shown 19 . In this regard, understanding how the effective energy of s-MTJs is characterised under STT, as well as magnetic field, is of significant interest not only from fundamental but also from technological aspects. The Néel-Arrhenius law under STT has been a longstanding question partly because the STT itself does not modulate the energy landscape due to its non-conservative nature, preventing one from defining its effective potential energy. Despite the difficulty, the expectation value of event time of the magnetisation switching, i.e. the Néel relaxation time τ, under field H and current I is phenomenologically expressed in a form: where τ 0 is~1 ns in magnetic systems 3 , and I C0 an intrinsic critical current. Regarding the exponent n I for the factor of current, different values, 1 (refs. 20,21 ) or 2 (refs. [22][23][24], have been theoretically derived, where the former was obtained by considering a fictitious temperature, whereas the latter was obtained by analysing the stochastic process based on the Fokker-Planck equation.
Experimentally, it has been practically inaccessible as far as one examines conventional thermally-stable MTJs and consequently, their decade-long unperturbed retention property has been extrapolated from limited data obtained in a reasonable time while assuming a certain number for n H or n I 25-33 . For applications with superparamagnetic tunnel junctions that actively utilise thermal fluctuation under STT, such uncertainty makes sophisticated engineering impractical as a rigorous description of modulation of the effective energy landscape is indispensable.
Here we experimentally study the Néel-Arrhenius law of a nanomagnet under STT utilising superparamagnetic tunnel junctions that allow direct determination of the event time under fields and currents [34][35][36][37][38][39] . Through measurements of homodynedetected ferromagnetic resonance (FMR) under current biases, high-speed STT switching with various pulse widths, and random telegraph noise (RTN) under various fields and currents, values of n H and n I are uniquely determined for given conditions. Furthermore, we show that by considering the local bifurcations under magnetic field and STT, which have not been considered for magnetic systems, the obtained results can be comprehended with the effects of the torques of the field and current without the difficulties to define the effective potential of STT.
Sample preparation and strategy of following experiments As shown in Fig. 1a, a stack structure, Ta(5)/Pt (5) 18 . Resistance (R)-area (A) product, RA, of the MgO tunnel barrier is 5.5 Ωμm 2 . The stack is patterned into MTJ devices, followed by annealing at 300°C. The MTJ device we will mainly focus on hereafter (device A) has a diameter D of 34 nm (results for device B with t CoFeB = 1.82 nm, RA = 8.1 Ωμm 2 , and D = 28 nm will be also shown later). In this size range, the magnetisation can be represented by a single vector without significant effects of spatial inhomogeneity for the present stack structure 16,27,40 . Both CoFeB layers have a perpendicular easy axis. Figure 1b shows the junction resistance R as a function of the perpendicular magnetic field H z . Gradual variation of mean R with H z and scattering of data points at the transition region reflects a superparamagnetic nature of the MTJ whose switching time is shorter than the measurement time of R (~1 s).

Ru(5) Ta(5) CoFeB
MgO To determine n H and n I in actual MTJs, we take into account the following two effects: electric-field modulation of magnetic anisotropy [41][42][43][44] and uncompensated stray field H S from the reference layer. Consequently, Δ, the argument of the exponential in Eq. (1), is rewritten as ð Þ ðVÞ V=V C0;PðAPÞ . V C0,P(AP) denotes the intrinsic critical voltage for STT switching from parallel, P, to antiparallel, AP, states (AP to P states). Because the electric-field effect on anisotropy is governed by the applied voltage V, we use an expression based on voltage input rather than current input. Equation (2) contains so many unknown variables (E 0 , H S , H K eff (V), V C0,P(AP) , n H and n I ) that one cannot directly determine the exponents by only measuring RTN. Thus, in the following, we first separately determine H K eff (V) from a homodyne-detected FMR and the next V C0 from the STT switching probability. After that, we determine H S and the exponents n H and n I from RTN measurement 34-38 as a function of V.
Electric-field effect on anisotropy field Firstly, we determine H K eff (V) from homodyne-detected FMR under dc bias voltage 45,46 . Figure 2a shows the circuit configuration for the measurement. Homodyne-detected voltage spectra are measured while sweeping H z at various frequencies and H K eff is determined from the peak position (see Methods and Supplementary Information for details). We perform this measurement at various dc biases and obtain H K eff vs. V, as shown in Fig. 2b. H K eff changes nonlinearly with V. We fit a quadratic equation to the obtained dependence and determine the coefficients for constant, linear, and quadratic terms to be μ 0 H K eff (0) = 77.0 ± 0.5 mT, μ 0 dH K eff /dV = −57.8 ± 1.6 mT V −1 , and μ 0 d 2 H K eff / dV 2 = −49.9 ± 7.5 mT V −2 (μ 0 is the permeability of vacuum). Note that the constant term represents the magnetic anisotropy field at zero bias whereas the linear and quadratic terms mainly originate from the electric-field modulation of anisotropy and an effect of Joule heating, respectively. The determined coefficients will be used in the analysis of RTN later.
Intrinsic critical voltage for STT switching Secondly, we determine V C0 for the STT switching. Because the fluctuation timescale of the studied system is on the order of milliseconds, we need to measure the junction state right after the switching pulse application when they are still non-volatile. To this end, we use a circuit configuration shown in Fig. 3a. This configuration is similar to that we used in our previous work to  Lines are linear fits, whose intercept yields intrinsic critical voltages V C0,P(AP) . study the switching error rate 47 ; a voltage waveform composed of initialisation, write, and read pulses [ Fig. 3b] is applied to the MTJ by an arbitrary waveform generator, and the transmitted signal at the read pulse is monitored by a high-speed oscilloscope to identify the final state of magnetisation configuration. The typical transmitted signals for P and AP states are shown in Fig. 3c. A clear difference is observed in the amplitude of the transmitted signal for different configurations due to the tunnel magnetoresistance. Write and read pulses are separated by 30 ns, which is much shorter than the shortest relaxation time shown later (~0.3 ms), ensuring a sufficiently low read-error rate (unintentional switching probability before the read pulse) <10 −4 (see Methods). The waveform is applied 200 times repeatedly and switching probability is evaluated. Figure 3d shows the write pulse voltage dependence of the switching probability with different write pulse durations t pulse . The switching voltage V C , defined as the voltage at 50% probability, is plotted as a function of the inverse of t pulse in Fig. 3e. In the precessional regime (typically t pulse ≲ several nanoseconds) where the switching/non-switching is determined by an amount of transferred angular momenta, V C is known to linearly depend on t pulse −1 and the intercept yields V C0 20,48 . From a linear fitting, V C0,P and V C0,AP are obtained as 313 ± 45 mV and −247 ± 42 mV, respectively.
Random telegraph noise measurement to determine the switching exponents With the results above, we are now ready to determine the switching exponents, n H and n I , from RTN measurement under various V and H z . Figures 4a, b show the circuit configurations for the measurement with V~0 and |V| ≥ 25 mV, respectively. For V~0, we apply a small direct current I = 200 nA and monitor R by an oscilloscope connected in parallel to the MTJ and probe the temporal magnetisation configuration. For |V| ≥ 25 mV, we monitor the divided voltage at reference resistor R r serially connected to the MTJ using an oscilloscope. Note that R r (= 470 Ω) is set to be much smaller than R to prevent a change in the electric-field effect on the magnetic anisotropy between P and AP states. Figure 4c shows typical results of RTN with various H z and the definition of the magnetisation switching event time t. Figure 4d shows the distribution of the number of unique t for μ 0 H z = −30.5 mT. As expected, the exponential distribution is confirmed, i.e., the number of events ∝ τ −1 exp(t/τ), indicating that the fluctuation is characterised by a Poisson process. From the fitting, expectation values of the event time for P and AP states, i.e., the relaxation time τ P and τ AP , are obtained as a function of H z as shown in Fig. 4e. Subsequently, Δ P(AP) can be determined from the relation τ = τ 0 expΔ P(AP) .
We measure τ P(AP) and Δ P(AP) for various H z and V. Figure 5a shows the obtained Δ P and Δ AP as a function of H z for various V. Δ P(AP) increases (decreases) with increasing H z for each V, as expected from the energy landscape modulation by H z . Also, the mean Δ gradually decreases with increasing V, which is also consistent with the trend of H K eff shown in Fig. 2b. To derive n H and n I , we then take the natural logarithm of the ratio between Δ P and Δ AP , which can be expressed from Eq. (2) as ln(Δ P /Δ AP ) vs. H z for each V is plotted in Fig. 5b. At a small perturbation limit, i.e., h, v P , v AP ≪ 1, Eq. (3) is reduced to ln(Δ P /Δ AP ) = 2n H h + n I (v AP − v P ); thus, the slope and intercept of the linear fit to the data shown in Fig. 5b give n H and n I , Fig. 4 Random telegraph noise. a, b Circuit configuration to measure random telegraph noise (RTN) with V~0 and V ≥ 25 mV, respectively. c Typical RTN signal monitored at the oscilloscope. d Histogram of event time (duration between two sequential switching events) for P and AP states for μ 0 H z = −30.5 mT. e Expected switching time as a function of perpendicular magnetic field H z . respectively. One can see that the results are well fitted by the linear function, validating the employed model. We analyse Fig. 5 with Eq. (3) and obtain n H and n I as a function of V for device A as shown in Fig. 6a. One can see that both n H and n I show similar values at each V and gradually decreases to about 1.5 with decreasing V. We perform the same procedure for the device B, whose properties are determined as μ 0 H K eff (0) = 129.0 ± 0.7 mT, μ 0 dH K eff /dV = −61.7 ± 2.3 mT V −1 , μ 0 d 2 H K eff /dV 2 = −58 ± 13 mT V −2 , V C0,P = 672 ± 4 mV and V C0,AP = −541 ± 2 mV. The obtained n H and n I are shown in Fig. 6b. At V = 0, the two devices show the same value for n H within experimental inaccuracy. Also, both n H and n I of device B show similar values with each other as in device A. However, in contrast to device A, they do not show meaningful variations at around 2 with V.

Discussion
As shown above, we have found that n H and n I show virtually the same value with each other for both devices A and B. Also, they are almost constant at around 2.0 for device B whereas change from 2.0 to 1.5 with V for device A. The main difference between devices A and B is t CoFeB , which manifests in a difference in μ 0 H K eff (0) [77.0 ± 0.5 mT for device A and 129.0 ± 0.7 mT for device B]. In the following, we will discuss the mechanism that can account for the obtained results in the context of the energy landscape and its bifurcation.
In systems with uniaxial anisotropy where the magnetic field is applied along the easy axis where the macrospin approximation holds, Brown derived n H = 2 in a high barrier region, E 0 ≳ k B T, using Kramers' analysis on the Fokker-Planck equation that is equivalent to the Landau-Lifshitz-Gilbert (LLG) equation with the Langevin term 3 . Taking into consideration the second-order anisotropy, where the magnetic anisotropy energy density is given by E = K 1 eff sin 2 θ + K 2 sin 4 θ, n H was pointed out to vary with K 2 / K 1 eff 7 , where K 1 eff , K 2 and θ are the firstand second-order effective anisotropy fields and polar angle of magnetisation vector, respectively. In the CoFeB/MgO system, positive voltage, which decreases electron density at the interface, was found to increase K 1 eff , while keeping μ 0 H K2 (≡μ 0 K 2 /4M S ) constant at around 45 mT 46,49,50 . Accordingly, as shown in the upper axes of Fig. 6a, b, in the present cases, K 2 /K 1 eff is calculated to be around 0.22 for device B whereas it increases up to 0.45 for device A. The numerical calculation, assuming material parameters of magnetic recording media (Δ 0 ≳ 60), shows that n H decreases from 2.0 to 1.5 in the range of K 2 /K 1 eff from 0 to 0.25 7 . In general, a dynamical system with pitchfork bifurcation leads to the switching exponent of n x = 2, while saddle-node bifurcation results in n x = 3/2. Note that the aforementioned magnetic energy density E gives the LLG equation dθ/dt = f(θ,H z ) = −αγμ 0 [(2K 1 eff /M S ) cosθ − (4K 2 /M S )cos 3 θ − H z ]sinθ. We show f(θ,H z ) takes two types of local bifurcations: pitchfork bifurcation appears at K 2 / K 1 eff < 0.25, while saddle-node bifurcation at K 2 /K 1 eff > 0.25 as shown in Fig. 6c, d, respectively [see Supplementary Information in detail]. Thus, the experimentally observed transition of the switching exponents is attributed to the transition of the bifurcation of the potential landscape through the modulation of K 2 / K 1 eff . However, the experiment shows the transition of n H at K 2 / K 1 eff ≈ 0.45, which is larger than that expected by the macrospin model (K 2 /K 1 eff = 0.25). This deviation implies that the local bifurcation of the magnetic potential and the resultant n H in the real MTJ device is more insensitive to higher-order anisotropy field K 2 than the macrospin limit, for example, due to the micromagnetic effects.
Regarding n I , some theoretical studies derived 1 by considering a fictitious temperature in LLG equation with the Langevin term 20,21 , whereas others derived 2 from an analysis of the Fokker-Planck equation [22][23][24] . Matsumoto et al. pointed out that n I rapidly decreases from 2 to 1.4 with increasing K 2 /K 1 eff from 0 to~0.25 51 . Experimentally, some assumed 1 25,[27][28][29][30]36,39 whereas others assumed 2 26,[31][32][33] , and importantly no studies access the number. The present experimental results support the scenario of Matsumoto et al., but, similarly to n H , the reduction of n I is more moderate than the theoretical prediction. This fact indicates that the mechanism for n H could be also applicable for the case with STT perturbation as well. Another important implication of our results is that, despite the non-conservative nature of STT, the pseudo energy landscape under STT can be investigated through the switching exponents. LLG equation with STT τ STT can be represented as dθ/dt = f(θ,x) = {−αγμ 0 [(2K 1 eff /M S )cosθ − (4K 2 / M S )cos 3 θ] + τ STT }sinθ. Since n H and n I show virtually the same value for all K 2 /K 1 eff conditions, meaning that f(θ,τ STT ) takes the same local bifurcation type as that for f(θ,H), our experiment reveals that in MTJ devices with perpendicular easy axis, the magnetic field and the STT effectively similarly modulate the energy landscape.
In summary, this work has experimentally revealed the hitherto-inaccessible representation of thermally-activated switching rate under field and STT, using a relevant material system for applications. The obtained results could allow for sophisticated engineering of non-volatile memory and unconventional computing hardware. Through the switching exponents, we have also accessed the local bifurcation of energy landscape under STT, and have found that, despite the qualitative difference between magnetic field and STT, their effect on the energy landscape is equivalent in the case of perpendicular MTJ. This work has also demonstrated that superparamagnetic tunnel junctions and analysis of their local bifurcation can serve as a versatile tool to investigate unexplored physics relating to thermally-activated phenomena in general with various configurations and external perturbations.  size determined from scanning electron microscopy observation and measured resistance for large devices with diameter D > 45 nm. The resistance-area product of device A (device B) is 5.5 Ωμm 2 (8.1 Ωμm 2 ), and the tunnel magnetoresistance ratio is 73% (74%). The nominal thickness of the MgO is 1.0 nm for both devices, and the difference of the RA corresponds to the~7% variation of the actual thickness due to the process variations between the two runs. D of devices A and B are determined from their resistance and RA to be D = 34 and 28 nm, respectively.

Methods
Homodyne-detected ferromagnetic resonance (FMR). With the circuit shown in Fig. 2a, homodyne-detected voltage-H z spectra were measured at various frequencies. As shown in previous papers, the spectra were well fitted by the Lorentz function and peak position was determined by the fitting 45,46 . From resonance frequency f r vs. H z , the effective anisotropy field H K eff was determined while assuming a constant second-order anisotropy field μ 0 H K2 = 45 mT 46,49,50 (μ 0 is the permeability of vacuum). The measurement was performed at various dc biases at AP configuration and obtained H K eff vs. V, as shown in Fig. 2b. A quadratic equation was fitted to the obtained dependence and the coefficients for constant, linear and quadratic terms were determined. Note that, in the error of H K eff , we have included the effect of H K eff difference in P and AP configurations due to the different device resistances and resultant Joule heatings in these configurations under the identical bias voltage.
Switching probability measurement. With the circuit shown in Fig. 3a, the switching probability was measured as functions of write pulse voltage amplitude and duration to determine intrinsic critical voltage V C0 . A voltage waveform composed of initialisation (0.45 V/300 ns), write (amplitude V write /duration t pulse ), and read (V read = 0.15 V/75 ns) pulses as shown in Fig. 3b, was generated by an arbitrary waveform generator (AWG). Both the interval of initialisation/write and write/read pulses were 30 ns, which is much shorter than the shortest relaxation time measured here (~0.3 ms), ensuring a sufficiently low read-error rate due to unintentional switching probability before the read pulse, exp(−30 ns/0.3 ms) ≤ 10 −4 . Single-shottransmitted voltage for write pulse was monitored to determine the magnetisation configuration; the transmitted voltage is~2Z 0 V read /(R + Z 0 ), where Z 0 is characteristic impedance 50 Ω, and due to the tunnel magnetoresistance, the transmitted voltage changes with magnetisation configuration. The typical transmitted signal for P and AP states is shown in Fig. 3c. Transmitted signals for 5 ns (between 15 and 20 ns in Fig. 3c) were averaged. Its averaged value 〈V〉 and standard error 〈(V-〈V〉) 2 〉 0.5 /N 0.5 for P and AP states were 2.44 ± 0.03 mV and 1.40 ± 0.03 mV, respectively (N is averaged points; 20 Gbit/s × 5 ns duration = 100 points), ensuring low read-error rate due to misassignment of the magnetisation configuration 47 . As shown in Fig. 3d, switching probability as a function of the voltage amplitude V at MTJ with write pulse duration t pulse from 1 to 5 ns was measured. The probability of the switching was determined from 200 times measurement. The switching measurement was conducted under H z of the stray field H S which was determined from the random telegraph noise measurement. Note that the anomaly of the switching probability at Ṽ −900 mV can be attributed to a change of the magnetic easy axis through the electric-field effect on the magnetic anisotropy, which is reported in previous works 52 . In addition, the slope of the switching probability at P sw~0 .5 for P to AP switching increases with decreasing the pulse duration, which is opposite to the thermally-stable MTJs. The decreases of the effective field and the thermal stability factor through the electric-field effect on magnetic anisotropy reasonably explain the behaviour.
Random telegraph noise (RTN). With the circuit shown in Fig. 4a, the RTN signal of the MTJs for V~0 was measured. Small direct current I = 200 nA was applied and R was monitored by an oscilloscope connected in parallel to the MTJ to probe the temporal magnetisation configuration. The voltage applied to MTJ here was up to 2.5 mV Fig. 6 Experimentally determined switching exponents n H n I . a, b Field-and current-induced switching exponents n H and n I as a function of bias voltage V for a device A and b device B. Error bars show standard error propagated from ferromagnetic resonance measurement and critical current measurement. The green band is a guide to the eye. c, d Bifurcation diagrams of the dθ/dt projected upon (θ,x)(= (θ,H) or (θ,I)) for c K 2 /K 1 eff < 0.25 and d K 2 /K 1 eff > 0.25, where K 1 eff , K 2 , θ, t, H, and I are the firstand second-order effective anisotropy fields, polar angle of magnetisation vector, the time, perpendicular magnetic field, and current. e, f Schematics of energy barrier of the effective magnetic potential corresponding to the yellow region in c and d, respectively.
(5 mV) for device A (device B), which is small enough to prevent major voltage/currentinduced effects in MTJs focused here. If one utilises the same circuit, the applied voltage for P and AP states varies by a factor of about 1.7 due to the tunnel magnetoresistance effect. Therefore, to prevent variation of effective anisotropy field for P and AP states, the circuit shown in Fig. 4b was utilised for |V| ≥ 25 mV. Direct voltage V to the MTJ was applied and divided voltage at reference resistor R r connected in serial to the MTJ was monitored using the oscilloscope. For measuring RTN on device A (device B), with setting R r = 0.47 kΩ (1 kΩ) much smaller than R, variation of applied voltages between P and AP states was prevented.
Attempt frequency. In determining Δ with the random telegraph noise measurement, Néel-Arrhenius raw τ P(AP) = τ 0 expΔ P(AP) with attempt frequency τ 0 of 1 ns was assumed. This assumption is widely adopted because τ is an exponential function of Δ and the value of τ 0 does not affect the estimated Δ, and τ 0 ranges between 0.1 and 10 ns. According to Brown's calculation with Kramer's method on the Fokker-Planck equation 3 , the attempt frequency τ 0 of the magnetic materials with uniaxial perpendicular magnetic anisotropy is [2αγμ 0 H K eff (1 − h 2 )(1 + h)] −1 (π/Δ 0 ) 0.5 under large barrier approximation Δ 0 )1, where α, γ, and μ 0 are damping constant 0.006, gyromagnetic ratio, and permeability of vacuum, respectively. In our devices, the Brown's attempt frequency above is derived to be 1.1 and 2.4 ns for device A and device B, respectively. Thus, the switching time τ vs. thermal stability factor Δ of device A should be well described by the Néel-Arrhenius law with τ 0 = 1 ns.

Data availability
The data that support the plots within this paper have been deposited in Zenodo at https://zenodo.org/record/6767828 53 .